//2D integration with y bound y1(x) and y2(x) and x bound x1 and x2  for func2d
#include "qgaus.h"
#include <math.h>

float cylr, sphr, dcsphacyl;
//---------------------------------------------------------------------------
struct NRf2   //For function of x and y
{
	double xsav;
	double (*func2d)(const double, const double);
	double operator()(const double y) //The integrand f(x;y) evaluated at fixed x
	{
	return func2d(xsav,y);//func2d(xsav,y);
	}
};
struct NRf1
{
	NRf2 f2;
	double (*y1)(double);
	double (*y2)(double);
	NRf1(double yy1(double), double yy2(double)) : y1(yy1),y2(yy2) {}
	double operator()(const double x) //This is H of eq. (4.8.5).
	{
	f2.xsav=x;
	return  qgaus(f2,y1(x),y2(x));

	}
};

template <class T>
double quad2d(T &func, const double x1, const double x2, double y1(double), double y2(double))
//Returns the integral of a user-supplied function func over a three-dimensional region specified
//by the limits x1, x2, and by the user-supplied functions y1, y2, z1, and z2, as defined in (4.8.2).
//Integration is performed by calling qgaus recursively.
{
NRf1 f1(y1,y2);
f1.f2.func2d=&func;
return qgaus(f1,x1,x2);
}
//The function to be integrated

//test function to be integrated 30xy with x=(0,1) and y=(x^2,x) ans=1.25
/*double func(const double x, const double y){return (30*x*y);}
//the equivalent functor, the functions defining the boundary :
double y1(const double x){return x*x;}
double y2(const double x){return x;}*/

/*
double func(const double x, const double y){return (x*x + y*y);}
//or as the equivalent functor, the functions defining the boundary can only be functions:
double y1(const double x){return -sqrt(1 - x*x);}
double y2(const double x){return sqrt(1 - x*x);}
//Integration of r^2 over a circle using quad2d_t
/*
	 Radius         Approx         Actual
   1.000000e-01   1.572260e-04   1.570796e-04
   2.000000e-01   2.515617e-03   2.513274e-03
   3.000000e-01   1.273531e-02   1.272345e-02
   4.000000e-01   4.024987e-02   4.021239e-02
   5.000000e-01   9.826628e-02   9.817477e-02
   6.000000e-01   2.037649e-01   2.035752e-01
   7.000000e-01   3.774997e-01   3.771482e-01
   8.000000e-01   6.439979e-01   6.433982e-01
   9.000000e-01   1.031560e+00   1.030599e+00
   1.000000e+00   1.572260e+00   1.570796e+00
   a=-1,b=1;
   ShowMessage("2D-INTEGRAL="+FloatToStr(quad2d(func, a, b,y1,y2)));
   */
double func(const double x, const double y)
{
		return sqrt(fabs(sphr*sphr - x*x - y*y));
}
//or as the equivalent functor, the functions defining the boundary can only be functions:
double lymin(const double x){return 0;}
double lymax(const double x){return GetMin(sqrt(fabs(cylr*cylr - (x-dcsphacyl)*(x-dcsphacyl))),sqrt(fabs(sphr*sphr - x*x)));}
